####################################
# Appendix H — Alternative RDD Windows
# Outputs: rdd_ideo_bin_year_rob1.png ... rdd_ideo_bin_year_rob5.png
# Panels:  (a) ±140, (b) ±130, (c) ±120, (d) ±110, (e) ±100
####################################

rm(list = ls())

suppressPackageStartupMessages({
  library(ggplot2)
  library(rdrobust)
  library(lfe)
  library(broom)
  library(rdd)   # for kernelwts()
})

# --- Directories ---
base_dir <- "~/Dropbox/Issue voting Chile"
fig_dir  <- file.path(base_dir, "09_replication/output")
data_dir <- file.path(base_dir, "09_replication/data")
dir.create(fig_dir, showWarnings = FALSE, recursive = TRUE)

# --- Mapping rob datasets -> panel/window labels & output paths ---
rob_ids  <- 1:5
panels   <- letters[1:5]  # (a)–(e)
windows  <- c("±140 days","±130 days","±120 days","±110 days","±100 days")

rob_rdata <- file.path(data_dir, sprintf("rdd_rob%d_2025aug19.RData", rob_ids))
out_png   <- file.path(fig_dir,  sprintf("rdd_ideo_bin_year_rob%d.png", rob_ids))

# (Optional) quick reference table
h_map <- data.frame(panel = panels, window = windows, rob = paste0("rob", rob_ids),
                    rdata = rob_rdata, png = out_png, stringsAsFactors = FALSE)
print(h_map)

# --- Helper to run one rob dataset and save the figure ---
run_one <- function(rdata_path, panel_label, window_label, outfile) {
  # Load objects: expected d06, d10, d14, d18 with vars:
  # ideo_bin, margin, treatment, cep, comuna
  env <- new.env(parent = emptyenv())
  load(rdata_path, envir = env)
  
  need <- c("d06","d10","d14","d18")
  if (!all(need %in% ls(env))) {
    stop(sprintf("Missing required data frames in %s", rdata_path))
  }
  
  get_piece <- function(df_name) {
    df <- get(df_name, envir = env)
    
    # Bandwidth selector
    bw <- rdrobust::rdbwselect(y = df$ideo_bin, x = df$margin)$bws[1]
    
    # Triangular kernel weights within bw
    w  <- rdd::kernelwts(df$margin, 0, bw = bw, kernel = "triangular")
    
    # FE RDD with interaction; FE: cep + comuna; cluster: comuna
    fit <- lfe::felm(
      ideo_bin ~ treatment + margin + treatment*margin | cep + comuna | 0 | comuna,
      data = df, weights = w, keepCX = TRUE
    )
    
    broom::tidy(fit, conf.int = TRUE, conf.level = 0.95)
  }
  
  r06 <- get_piece("d06")
  r10 <- get_piece("d10")
  r14 <- get_piece("d14")
  r18 <- get_piece("d18")
  
  # Extract treatment coefficient (first row, mirroring your script)
  pe      <- c(r06$estimate[1], r10$estimate[1], r14$estimate[1], r18$estimate[1])
  se      <- c(r06$std.error[1], r10$std.error[1], r14$std.error[1], r18$std.error[1])
  upper95 <- c(r06$conf.high[1], r10$conf.high[1], r14$conf.high[1], r18$conf.high[1])
  lower95 <- c(r06$conf.low[1],  r10$conf.low[1],  r14$conf.low[1],  r18$conf.low[1])
  
  years <- factor(c("1996-2006","1996-2010","1996-2014","1996-2018"),
                  levels = c("1996-2006","1996-2010","1996-2014","1996-2018"))
  results <- data.frame(years, pe, se, lower95, upper95)
  
  panel_title <- sprintf("(%s) %s", panel_label, window_label)
  
  p <- ggplot(results, aes(x = years, y = pe)) +
    geom_pointrange(aes(ymin = lower95, ymax = upper95)) +
    scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1)) +
    coord_cartesian(ylim = c(-0.2, 1)) +
    geom_hline(yintercept = 0) +
    xlab("Years") + ylab("RDD estimates") +
    ggtitle(panel_title) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5),
          legend.position = "none",
          text = element_text(size = 20))
  
  ggsave(filename = outfile, plot = p, width = 10, height = 7)
  message("Saved: ", outfile)
}

# Figure H1

# --- Produce all five panels ---
for (i in seq_along(rob_rdata)) {
  run_one(rob_rdata[i], panels[i], windows[i], out_png[i])
}
